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ABSTRACT 

We investigate the evolution of non-linear density perturbations by taking into ac- 
count the effects of deviations from spherical symmetry of a system. Starting from 
the standard spherical top hat model in which these effects are ignored, we introduce 
a physically motivated closure condition which specifies the dependence of the addi- 
tional terms on the density contrast, S. The modified equation can be used to model 
the behaviour of an overdense region over a sufficiently large range of 5. The key new 
idea is a Taylor series expansion in (1/5) to model the non-linear epoch. We show that 
the modified equations quite generically lead to the formation of stable structures in 
which the gravitational collapse is halted at around the virial radius. The analysis 
also allows us to connect up the behaviour of individual overdense regions with the 
non-linear scaling relations satisfied by the two point correlation function. 

Key words: Cosmology : theory - dark matter, large scale structure of the Universe 



1 INTRODUCTION 

Analytic modelling of the non-linear phase of gravitational 
clustering has been a challenging but interesting problem 
upon which a considerable amount of attention has been 
bestowed in recent years. The simplest, yet remarkably suc- 
cessful, model for non-linear evolution is the Spherical Col- 
lapse Model (SCM, hereafter), which has been applied in the 
study of various empirical results in the gravitational insta- 
bility paradigm. Unfortunately, this approach has serious 
flaws — both mathematically and conceptually. Mathemat- 
ically, the SCM has a singular behaviour at finite time and 
predicts infinite density contrasts for all collapsed objects. 
Conceptually, it is not advisable to model the real universe 
as a sphere, in spite of the standard temptations to which 
theoreticians often succumb. The two issues are, of course, 
quite related, since, in any realistic situation, it is the devi- 
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ations from spherical symmetry which lead to virialised sta- 
ble structures getting formed. In conventional approaches, 
this is achieved by an ad hoc method which involves halting 
the collapse at the virial radius by hand and mapping the 
resulting non-linear and linear overdensities to each other. 
This leads to the well known rule-of-thumb that, when the 
linear overdensity is about 1.68, bound structures with non- 
linear overdensities of about 178 would have formed. The 
singular behaviour, however, makes the actual trajectory of 
a spherical system quite useless after the turnaround phase 
— a price we pay for the arbitrary procedure used in sta- 
bilizing the system. But the truly surprising feature is that, 
despite its inherent arbitrariness, the SCM, when properly 
interpreted, seems to give useful insights into the behaviour 
of real systems. The Press-Schechter formalism (Press & 
Schechter 1974), for the abundance of bound structures, 
uses SCM implicitly; more recently, it was shown that the 
basic physics behind the non-linear scaling relations (NSR) 
obeyed by the two point correlation function can be obtained 
from a judicious application of SCM (Padmanabhan 1996a). 
These successes, as well as the inherent simplicity of the un- 
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derlying concepts, make the SCM an attractive paradigm 
for studying non-linear evolution in gravitational clustering 
and motivate one to ask : Can we improve the basic model 
in some manner so that the behaviour of the system after 
turnaround is 'more reasonable' ? 

It is clear from very general considerations that such an 
approach has to address fairly non-trivial technical issues. 
To begin with, exact modelling of deviations from spheri- 
cal symmetry is quite impossible since it essentially requires 
solving the full BBGKY hierarchy. Secondly, the concept of 
a radius R(t) for a shell, evolving only due to the gravita- 
tional force of the matter inside, becomes ill-defined when 
deviations from spherical symmetry are introduced. Finally, 
our real interest is in modelling the statistical features of the 
density growth; whatever modifications we make to SCM 
should eventually tie up with known results for the evolu- 
tion of, for instance, the two point correlation function. That 
is, we have to face the question of how best to obtain the 
statistical properties of the density field from the behaviour 
of a single system. 

In this paper, we try to address these problems in a limited 
but focussed manner. We tackle the deviations from spher- 
ical symmetry by the retention of a term (which is usually 
neglected) in the equation describing the growth of the den- 
sity contrast. Working in the fluid limit, we show that this 
term is physically motivated and present some arguments to 
derive an acceptable form for the same. The key new idea is 
to introduce a Taylor series expansion in (1/5) (where 5 is 
the density contrast) to model the non-linear evolution. We 
circumvent the question of defining the 'radius' of the non- 
spherical regions by working directly with density contrasts. 
Finally, we attempt to make the connection with statisti- 
cal descriptors of non-linear growth, by using the non-linear 
scaling relations known from previous work. More precisely, 
we show that the modified equations predict a behaviour 
for the relative pair velocity (when interpreted statistically) 
which agrees with the results of N-body simulations. 
The paper is divided into the following sections. The rele- 
vant equations describing the SCM are set out in Section 
2; we also summarise the physical and ad hoc aspects of 
the SCM here. Next, we recast the equations in a different 
form and introduce two functions (i) a "virialization term" 
and (ii) a function /isc(5), whose asymptotic forms are easy 
to determine. The behaviour of hsc{5) in the presence and 
absence of the "virialization term" is also detailed here. In 
Section 4, we present the arguments that give the functional 
forms for the above term over a large range of 5; we then go 
on to present the results in terms of a single collapsing body 
and show how this term stabilizes a collapse which would 
have otherwise ended up in a singularity in terms of the 
growth of the density contrast with time. When this term 
is carried through into the equation for R(t) for a single 
system, it can be seen the radius reaches a maximum and 
gracefully decreases to a constant, remaining so thereafter. 
In the standard SCM, the radius decreases from the maxi- 
mum all the way down to zero, thereby causing the density 



to diverge. Section 5 summarises the results and discusses 
their implications. 



2 THE SPHERICAL COLLAPSE MODEL. 

The scales of interest in the current work are much smaller 
than the Hubble length and the velocities in question are 
non-relativistic; Newtonian gravity can hence be used for 
the following analysis. We will consider the case of a dust- 
dominated, Q — 1 Universe and treat the system in the fluid 
limit as being made up of pressureless dust of dark mat- 
ter, with a smoothed density, p m (t,x.), and a mean veloc- 
ity, v(i,x). (This approach, of course, ignores effects arising 
from shell crossing and multi-streaming; these will be com- 
mented on later.) 

The density contrast, 5(x, i) is defined by 



p m (t,x) = Pb (t)[l + 6(x,t)] 



(1) 



where pt denotes the smooth background density of matter. 
We define a velocity field u 1 = v l /(ad), where v z is the pe- 
culiar velocity (obtained after subtracting out the Hubble 
expansion) and a(t) denotes the scale factor. Taking the di- 
vergence of the field u % and writing it as 



(2) 



where Oij is the shear tensor, Q is the rotation vector and 9 
is the expansion, we can manipulate the fluid equations (see 
e.g. Padmanabhan 1996b) to obtain the following equation 
for 5 



d 2 8 3 d8 3 ... 
d^ + Tad- a -2a-i 5{1 + 5) 



4 i . 2 

3(1 + 5) 

The same equation can be written in terms of time t as 



■■ 4 8 2 
o — — 



3(1 + 5) + a 



4nGp b 8(l + 8) +d 2 (l + 8)(a 2 - 2fl 2 ) (4) 



This equation turns out to be the same as the one for den- 
sity contrast in the SCM, except for the additional term 
(1 + S)(a 2 — 2f2 2 ), arising from the angular momentum and 
shear of the system. To see this explicitly, we introduce a 
function R(t) by the definition 



1 + 5 



9GMt 



X- 



(5) 



2R 3 R 3 

where M and A are constants. Using this relation between 5 
and R(t), equation (^) can be converted into the following 
equation for R(t) 
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k = - q tf-¥ (° 2 -™ 2 )R (6) 

where the first term represents the gravitational attraction 
due to the mass inside a sphere of radius R and the second 
gives the effect of the shear and angular momentum. 
In the case of spherically symmetric evolution, the shear and 
angular momentum terms can be set to zero; this gives 

d 2 R _ _GM_ 

dt 2 ~ R 2 V> 

which governs the evolution of a spherical shell of radius R, 
collapsing under its own gravity; M can now be identified 
with the mass contained in the shell; this is standard SCM. 

At this point, it is important to note a somewhat subtle 
aspect of these equations. The original fluid equations are 
clearly Eulerian in nature: i.e. the time derivatives give the 
temporal variation of the quantities at a fixed point in space. 
However, the time derivatives in equation (^), for the den- 
sity contrast 8, are of a different kind. Here, the observer 
is moving with the fluid element and hence, in this, La- 
grangian case, the variation in density contrast seen by the 
observer has, along with the intrinsic time variation, a com- 
ponent which arises as a consequence of his being at different 
locations in space at different instants of time. When the 5 
equation is converted into an equation for the function R(t), 
the Lagrangian picture is retained; in SCM, we can interpret 
R(t) as the radius of a spherical shell, co-moving with the 
observer. The mass M within each shell remains constant in 
the absence of shell crossing (which does not occur in the 
standard SCM for reasonable initial conditions) and the en- 
tire formalism is well defined. The physical identification of 
R is, however, not so clear in the case where the shear and 
rotation terms are retained, as these terms break the spher- 
ical symmetry of the system. We will nevertheless continue 
to think of R as the "effective shell radius" in this situation, 
defined by equation (^) governing its evolution. Of course, 
there is no such ambiguity in the mathematical definition of 
R in this formalism. 

Before proceeding further, let us briefly summarize the re- 
sults of standard SCM. Equation (R) can be integrated to 
obtain R(t) in the parametric form 

R = cos 0) (8) 

where Ri, Si and t% are the initial radius, initial den- 
sity contrast and initial time, respectively, with R\ = 
(9GMt?/2)(l + Si)' 1 ~ (9GMt?/2) for S t < 1. Given M, 
there are only two independent constants, viz U and Si. 
All the physical features of the SCM can be easily derived 
from the above solution. Each spherical shell expands at a 
progressively slower rate against the self-gravity of the sys- 



tem, reaches a maximum radius and then collapses under 
its own gravity, with a steadily increasing density contrast. 
The maximum radius, Rmax = Ri/Si, achieved by the shell, 
occurs at a density contrast 5 = (9n 2 /16) - lfs 4.6, which 
is in the "quasi-linear" regime. In the case of a perfectly 
spherical system, there exists no mechanism to halt the in- 
fall, which proceeds inexorably towards a singularity, with 
all the mass of the system collapsing to a single point. Thus, 
the fate of the shell (as described by equations (^) and is 
to collapse to zero radius at 6 — 2tt with an infinite density 
contrast; this is, of course, physically unacceptable. 
In real systems, however, the implicit assumptions that (i) 
matter is distributed in spherical shells and (ii) the non- 
radial components of the velocities of the particles are small, 
will break down long before infinite densities are reached. In- 
stead, we expect the collisionless dark matter to reach virial 
equilibrium. After virialization, \U\ = 2K, where U and K 
are, respectively, the potential and kinetic energies; the virial 
radius can be easily computed to be half the maximum ra- 
dius reached by the system. 

The virialization argument is clearly physically well- 
motivated for real systems. However, as mentioned earlier, 
there exists no mechanism in the standard SCM to bring 
about this virialization; hence, one has to introduce by hand 
the assumption that, as the shell collapses and reaches a par- 
ticular radius, say Rmax/2, the collapse is halted and the 
shell remains at this radius thereafter. This arbitrary intro- 
duction of virialization is clearly one of the major drawbacks 
of the standard SCM and takes away its predictive power in 
the later stages of evolution. We shall now see how the re- 
tention of the angular momentum term in equation (^|) can 
serve to stabilize the collapse of the system, thereby allowing 
us to model the evolution towards r V i r = R ma x/2 smoothly. 



3 THE Hsc(S) FUNCTION. 

As detailed in the previous section, the primary defect of 
the standard SCM is the ad hoc nature of the stabilization 
of the shell against its collapse under gravity, which arises 
on account of the assumption of perfect spherical symmetry, 
implicit in the neglect of the shear and angular momentum 
terms. We hence return to equation (^), retain the above 
terms, and recast the equation into a form more suitable 
for analysis. Using logarithmic variables, Dsc = In (1 + S) 
and a = In a, equation (||) can be written in the form (the 
subscript 'SC stands for 'Spherical Collapse') 

ri 2 £>sc 1 ( dP S c \ 2 1 dD sc 
da 2 3 V da J 2 da 

| [expose) - f] + a {a 2 - 2fi 2 ) (10) 
It is convenient to introduce the quantity, S, defined by 

S = a 2 {a 2 -2Q 2 ) (11) 
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which we shall hereafter call the "virialization term". The 
consequences of the retention of the virialization term are 
easy to describe qualitatively. We expect the evolution of 
an initially spherical shell to proceed along the lines of the 
standard SCM in the initial stages, when any deviations 
from spherical symmetry, present in the initial conditions, 
are small. However, once the maximum radius is reached and 
the shell recollapses, these small deviations are amplified by 
a positive feedback mechanism. To understand this, we note 
that all particles in a given spherical shell are equivalent due 
to the spherical symmetry of the system. This implies that 
the motion of any particle, in a specific shell, can be con- 
sidered representative of the motion of the shell as a whole. 
Hence, the behaviour of the shell radius can be understood 
by an analysis of the motion of a single particle. The equa- 
tion of motion of a particle in an expanding universe can be 
written as 



Xi + 2 — Xi — r- 

a a A 



(12) 



where a(t) is the expansion factor of the locally overdense 
"universe" . The Xi term acts as a damping force when it is 
positive; i.e. while the background is expanding. However, 
when the overdense region reaches the point of maximum 
expansion and turns around, this term becomes negative, 
acting like a negative damping term, thereby amplifying any 
deviations from spherical symmetry which might have been 
initially present. Non-radial components of velocities build 
up, leading to a randomization of velocities which finally re- 
sults in a virialised structure, with the mean relative velocity 
between any two particles balanced by the Hubble flow. It 
must be kept in mind, however, that the introduction of 
the virialization term changes the behaviour of the solution 
in a global sense and it is not strictly correct to say that 
this term starts to play a role only after recollapse, with the 
evolution proceeding along the lines of the standard SCM 
until then. It is nevertheless reasonable to expect that, at 
early times when the term is small, the system will evolve 
as standard SCM to reach a maximum radius, but will fall 
back smoothly to a constant size later on. 
The virialization term, S, is, in general, a function of a and 
x, especially since the derivatives in equation (|l|) are total 
time derivatives, which, for an expanding Universe, contain 
partial derivatives with respect to both x and t separately. 
Handling this equation exactly will take us back to the full 
non-linear equations for the fluid and, of course, no progress 
can be made. Instead, we will make the ansatz that the 
virialization term depends on t and x only through <5(i,x). 



S(a,x) = S((5(a,x)) = S(Z>sc) 



(13) 



In other words, S is a function of the density contrast alone. 
This ansatz seems well motivated because the density con- 
trast, 8, can be used to characterize the SCM at any point 
in its evolution and one might expect the virialization term 
to be a function only of the system's state, at least to the 
lowest order. Further, the results obtained with this assump- 



tion appear to be sensible and may be treated as a test of 
the ansatz in its own framework. 

To proceed further systematically, we define a function hsc 
by the relation 



dDsc 
da 



= 3ftsc 



(14) 



For consistency, we shall assume the ansatz ftsc(a,x) = 
hsc \S(a, x)]. The definition of ftsc allows us to write equa- 
tion (hot) as 



dhsc _ ,2 



hsC . 1 r m \ 11 , S(Dsc) ,.r\ 

— + 2 [exp(£>sc) - 1] + g (15) 



Dividing ( |l5| ) by (|l4|), we obtain the following equation for 
the function hsc(Dsc) 



dhs 



dDsc 



hsc 
3 



1 1 

6 6ftsc 



exp(D s 



H + 



S(D S 



9h sc 



(16) 



If we know the form of either hsc{Dsc) or S(Dsc), this 
equation allows us to determine the other. Then, using equa- 
tion (|l4j), one can determine Dsc- Thus, our modification 
of the standard SCM essentially involves providing the form 
of Ssc(Dsc) or ftsc(fsc)- We shall now discuss several fea- 
tures of such a modelling in order to arrive at a suitable 
form. 

The behaviour of hsc{Dsc) can be qualitatively understood 
from our knowledge of the behaviour of 8 with time. In the 
linear regime (8 <C 1), we know that <5 grows linearly with 
a; hence hsc increases with Dsc- At the extreme non-linear 
end (i5 > 1), the system "virializes" , i.e. the proper radius 
and the density of the system become constant. On the other 
hand, the density pb, of the background, falls like t~ 2 (or 
a~ s ) in a flat, dust-dominated Universe. The density con- 
trast is defined by <5 = (p/pb — 1) ~ p/pb (for i5 > 1) and 
hence 



8 oct oc a 



(17) 



in the non-linear limit. Equation (JIJ) then implies that 
hsc{8) tends to unity for i > 1. Thus, we expect that 
hsc(Dsc) will start with a value far less than unity, grow, 
reach a maximum a little greater than one and then 
smoothly fall back to unity. [A more general situation dis- 
cussed in the literature corresponds to ft — ► constant as 
8 — > oo, though the asymptotic value of ft is not necessarily 
unity. Our discussion can be generalised to this case and we 
plan to explore this in a future work.] 

This behaviour of the ftsc function can be given another 
useful interpretation whenever the density contrast has a 
monotonically decreasing relationship with the scale, x, with 
small x implying large 8 and vice-versa. Then, if we use a 
local power law approximation 8 oc x~ n for i > 1 with some 
n > 0, -Dsc oc ln(a; _1 ) and 
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, dDsc 
hsc oc — - — oc 
da 



din a 



xa 
oc — oc 
ax 



v 

ax 



(18) 



where v = ax denotes the mean relative velocity. Thus, hsc 
is proportional to the ratio of the peculiar velocity to the 
Hubble velocity. We know that this ratio is small in the lin- 
ear regime (where the Hubble flow is dominant) and later 
increases, reaches a maximum and finally falls back to unity 
with the formation of a stable structure; this is another ar- 
gument leading to the same qualitative behaviour of the hsc 
function. 

Note that, in standard SCM (for which S = 0), equation 
( [HI ) reduces to 



„, dhsc ,2 h sc . 5 
3hsC dD^ =hsc - — + 2 



(19) 



The presence of the linear term in 8 on the RHS of the above 
equation causes hsc to increase with 5, with hsc oc S 1 ^ 2 for 
8 S> 1. If virialization is imposed as an ad hoc condition, then 
hsc should fall back to unity discontinuously — which is 
clearly unphysical; the form of S(S) must hence be chosen so 
as to ensure a smooth transition in hsc(S) from one regime 
to another. 

As an aside, we would like to make some remarks on the na- 
ture of the "virialization term", S(5), in a somewhat wider 
context. As is well-known, gravitational clustering can be 
described at three different levels of approximation, by dif- 
ferent mathematical techniques. The first approach tracks 
the clustering by following the true particle trajectories; this 
is what is done, for example, in N-body simultations. This 
method does not involve any approximation (other than the 
validity of the Newtonian description at the scales of inter- 
est); it is, however, clearly analytically intractable. At the 
next level, one may describe the system by an one-particle 
distribution funtion and attempt to solve the collisionless 
Boltzmann equation for the distribution function /(t,x,v); 
the approximation here lies in the neglect of gravitational 
collisions, which seems quite reasonable as the time scale 
for such collisions is very large for standard dark matter 
particles. Finally, one can treat the system in the fluid limit 
described by five functions: the density p(t,x), mean veloc- 
ity v(t,x), and gravitational potential <j>(t,x), thus neglect- 
ing multi-streaming effects. Our analysis was based on this 
level of approximation. The key difference between the last 
two levels of description lies in the fact that the distribu- 
tion function allows for the possibility of different particle 
velocities at any point in space (i.e. the existence of velocity 
dispersions) , while the fluid picture assumes a mean velocity 
at each point. It is also known that the gradients in veloc- 
ity dispersion can provide a kinetic pressure which will also 
provide support against gravitational collapse. While a de- 
tailed analysis of these terms is again exceedingly difficult, 
one can incorporate the lowest order effects of the gradient 
in the velocity dispersion by modifying equation ( |l0| ) to the 
form 



d 2 D SC 1 ( dD S c \ 2 1 dP sc 
da 2 3 V da J 2 da 

| [expose) - 1] + a 2 {o 2 - 2Q 2 ) + f(a, x) (20) 

where f(a, x) contains the lowest order contributions from 
the dispersion terms. We can then define 



S(a,x) = a 2 (a 2 - 2fi 2 ) + f(a,x) 



(21) 



and again invoke the ansatz S(a,x) = S(5). Note that S(5) 
now contains the lowest order contributions arising from 
shell crossing, multi-streaming, etc., besides the shear and 
angular momentum terms, i.e. it contains all effects leading 
to virialization of the system. We demonstrate explicitly in 
an appendix that velocity dispersion terms arise naturally in 
the "force" equation (for the function h = —v/ax), derived 
from the BBGKY hierarchy, and play the same role as the 
function S(S) in the fluid picture. This clearly justifies the 
above procedure and shows that our approach could have a 
somewhat larger domain of validity than might be expected 
from an analysis based on the fluid picture. 



4 THE VIRIALIZATION TERM 

We will now derive an approximate functional form for the 
virialization function from physically well-motivated argu- 
ments. If the virialization term is retained in equation (pf), 
we have 



H 2 R 



(22) 



d 2 R _ GM_ 

~dt T ~ ~"R2 3 

where H = a/a. Let us first consider the late time behaviour 
of the system. When virialization occurs, it seems reasonable 
to assume that R constant and R — > 0. This implies that, 
for large density contrasts, 



3GM 



(S » 1) 



R 3 H 2 

Using H — a/a — (2/3i), and equation (|^) 

27GMt 2 3., „ 3 r 1N 



(23) 



(24) 



Thus, the "virialization" term tends to a value of (—35/2) in 
the non-linear regime, when stable structures have formed. 
This asymptotic form for S(5) is, however, insufficient to 
model its behaviour over the larger range of density con- 
trast (especially the quasi-linear regime) which is of interest 
to us. Since S(S) tends to the above asymptotic form at 
late times, the residual part, i.e. the part that remains af- 
ter the asymptotic value has been subtracted away, can be 
expanded in a Taylor series in (1/5) without any loss of gen- 
erality. Retaining the first two terms of expansion, we write 
the complete virialization term as 
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(25) 



Replacing for S(S) in equation (0), we obtain, for i > 1 
d/isc .2 , h S c 1 A B 

[It can be easily demonstrated that the first order term in the 
Taylor series is alone insufficient to model the turnaround 
behaviour of the h function. We will hence include the next 
higher order term and use the form in equation (^) for the 
virialization term. The signs are chosen for future conve- 
nience, since it will turn out that both A and B are greater 
than zero.] In fact, for sufficiently large 8, the evolution de- 
pends only on the combination q = (B/A 2 ). This is most 
easily seen by rewriting equation (^|), replacing S(S) with 
the above form. Taking the limit of large 8, i.e. 8 3> 1, and 
rescaling 8 to 8/ A, we obtain 



d 2 8 3 d8 



db 2 + 2b db 



4 

37 



(-) 

KdbJ 



a 
1 

—a + 
a 2 - 



(27) 
(28) 



B 1 
~A 2 "a T 8 

q 

a 2 8 

From the form of the equation it is clear that the constants 
A and B occur in the combination q — B /A 2 and hence the 
non-linear regime is modelled by a one parameter family for 
the virialization term. 
Equation (Eih can be written as 



7? = 



GM 



4R 
27*2 



27GMt 
4i? 3 



A 
J 



B 



(29) 



Using 8 = 9GMt 2 /2R 3 and B = qA 2 we may express equa- 
tion (^) completely in terms of R and t. We now rescale R 
and t in the form R = r v i r y(x) and t = fix, where r v i r is 
the final virialised radius [i.e. R —* r V i r for t oo], where 
(i 2 = (8/3 5 )(A/GM)rl ir , to obtain the following equation 
for y(x) 



y 



y 



27 



V 



4 9 i 6 



(30) 



We can integrate this equation to find a form for y q (x) 
(where y q (x) is the function y(x) for a specific value of q) 
using the physically motivated boundary conditions y — 1 
and y' = as x — > oo, which is simply an expression of 
the fact that the system reaches the virial radius r v i r and 
remains here thereafter. 

The results of numerical integration of this equation for a 
range of q values are shown in figure (Q). As expected on 
physical grounds, the function has a maximum and grace- 
fully decreases to unity for large values of x [the behaviour 
of y(x) near x — is irrelevant since the original equation 
is valid only for 8 > I , at least]. For a given value of q, it is 
possible to find the value x c at which the function reaches its 
maximum, as well as the ratio y m ax = Rmax/r V i r . The time, 
tmax, at which the system will reach the maximum radius is 
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Figure 1. The figure shows the function y q (x) for some values 
of q. The x axis has scaled time, x and the y axis is the scaled 
radius y. 



related to x c by the relation t max = fix c = <o(l + z max )~ i/2 , 
where to = 2/(3ifo) is the present age of the universe and 
Zmax is the redshift at which the system turns around. Fig- 
ure (Q) shows the variation of x c and y max = (R m ax/r V i r ) 
for different values of q. The entire evolution of the system 
in the modified spherical collapse model (MSCM) can be 
expressed in terms of 



R(t) = r viT y q {t/0) 



(31) 



where fi = (t /x c )(l + z max ) 3/2 . 

In SCM, the conventional value used for (r V i r /R„ lax ) is 
(1/2), which is obtained by enforcing the virial condition 
that \U\ = 2K, where U is the gravitational potential en- 
ergy and K is the kinetic energy. It must be kept in mind, 
however, that the ratio (r V i r / Rmax) is not really constrained 
to be precisely (1/2) since the actual value will depend on the 
final density profile and the precise definitions used for these 
radii. While we expect it to be around 0.5, some amount of 
variation, say between 0.25 and 0.75, cannot be ruled out 
theoretically. 

Figure (^) shows the parameter (Rmax/r V i r ), plotted as a 
function of q = B/A 2 (dashed line), obtained by numerical 
integration of equation ( |22[ ) with the ansatz (^5[). The solid 
line gives the dependence of x c (or equivalently t max ) on 
the value of q. It can be seen that one can obtain a suitable 
value for the (r V i r / R max ) ratio by choosing a suitable value 
for q and vice versa. Using equation ([lij) and the definition 
8 oc t 2 /R 3 , we obtain 



Aao(x) = l-|^ 
2 y dx 



(32) 
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Density contrast, <5 



Figure 2. The figure shows the parameters (Rmax /r v i r ) (broken 
line) and x c (solid line) as a function of q = B/A 2 . This clearly 
demonstrates that the single parameter description of the viri- 
alization term is constrained by the value that is chosen for the 

ratio V v i r J Rmax • 



which gives the form of hsc{x) for a given value of q; this, 
in turn, determines the function y q (x). Since S can be ex- 
pressed in terms of x, y and x c i is 8 = (9-K 2 /2x 2 c )x 2 /y 3 , this 
allows us to implicitly obtain a form for hsc{8), determined 
only by the value of q. Figure (Eh shows the behaviour of 
hsc functions obtained by integrating equation (jl^) back- 
wards, assuming that hsc — > 1 as S — » oo. It is seen that all 
the curves have the same turnaround behaviour expected on 
the basis of the physical arguments presented in the earlier 
section. 

If the functional form for hsc - determined, say, from N- 
body simulations - is used as a further constraint, we should 
be able to obtain the values of q. The major hurdle in at- 
tempting to do this is the fact that the available simulation 
results are given in terms of the averaged two point corre- 
lation function, £, and the averaged pair velocity, h(a,x), 
defined by 



£(x, a)x dx ; h(a, x) 



(v(a,x)) 



(33) 



where the two-point correlation function £ is defined as the 
Fourier transform of the power spectrum, P(k), of the distri- 
bution. The results published in the literature assume that 
h(a,x) depends on a and x only through £(a,x), that is, 
h(a,x) = h[£(a,x)]. This assumption has been invoked in 
several papers in the past (See e.g. Hamilton et al. 1991, 
Nityananda & Padmanabhan 1994, Mo et al. 1995, Padman- 
abhan 1996a, Padmanabhan & Engineer 1998) and seems to 
be validated by numerical simulations. The fitting formula 



Figure 3. The figure shows the /igc function, obtained for various 
values of q. The values of q and y m ax = Rmax/^vir for the curves 
are indicated at the top right hand corner. (Further discussion in 
text) 



for h(£) can be obtained from related fitting formulas avail- 
able in the literature (e.g. Hamilton et al. 1991). These are, 
however, statistical quantities and are not well defined for 
an isolated overdense region. Hence we have to first make 
the correspondence between hsc{8) and h(£), which we do 
as follows. 

It is possible to show by standard arguments (Nityananda 
& Padmanabhan 1994) that 



1=3,(1+1) 
that is, 



da 



(34) 



(35) 



where D — ln(l + £) and a = In a. Equation (|35|) is very 
similar to equation (|l4|) , which defines the function hsc(8), 
except for the different definitions of D and Dsc in terms 
of £ and S respectively. This suggests that one can obtain 
a relation between hsc(S) and /i(£) by relating the density 
contrast <5 of an isolated spherical region to the two-point 
correlation function £ averaged over the distribution at the 
same scale. We essentially need to find a mapping between 
£ and S which is valid in a statistical sense. 
Gravitational clustering is known to have three regimes in 
its growing phase, usually called "linear" , "quasi-linear" and 
"non-linear" respectively. The three regimes may be charac- 
terized by values of density contrast as 8 <C 1 in the linear 
regime, 1 < S < 100 in the quasi-linear regime and 100 < 5 
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in the non-linear regime. The three regimes have different 
rates of growth for various quantities of interest such as 8, £ 
and so on. In the linear regime, it is well known that the den- 
sity contrast grows proportional to the scale factor, a. This 
implies that the power spectrum, P(k) = \8k\ 2 (where 8k is 
the Fourier mode corresponding to 5(x)), grows as a 2 . Con- 
sequently, £, which is related to P(k) via a Fourier transform, 
also grows as a 2 , i.e. as the square of the density contrast. In 
the quasi-linear and non-linear regimes, the density contrast 
does not grow linearly with the scale factor and the relation 
between 8 and £ is not so clearly defined. The quasi- linear 
regime may be loosely construed as the interval of time dur- 
ing which the high peaks of the initial Gaussian random field 
have collapsed, although mergers of structures have not yet 
begun to play an important role. (This idea was used in Pad- 
manabhan (1996a) to model the non- linear scaling relations 
successfully). If we consider a length scale smaller than the 
size of the collapsed objects, the dominant contribution to 
£ (at this scale) arises from the density profiles centered on 
the collapsed peaks. Using the relation 



Pb 



(36) 



for density profiles around high peaks, one can see that f oc 8 
in this regime. In the non-linear regime, 8 and £ have the 
forms 5(a,x) — a [i F(ax), £ = a 3 G(ax), where i is a co- 
moving and r — ax is a proper coordinate. When the system 
is described by Lagrangian coordinates (which correspond to 
proper coordinates r = ax, i.e. at constant r), £ is propor- 
tional to 8. Thus, the relation £ oc 8 appears to be satisfied 
in all regimes, except at the very linear end. Since we are 
only interested in the 8 > 1 range, we use £ ~ <5 and compare 
equations (Q) and (^) to identify 



hsc{5) « h[i) 



(37) 



It is now straightforward to choose the value of q such that 
the known fitting function for the h function is reproduced 
as closely as possible. We use the original function given by 
Hamilton et al. (1991) to obtain the following expression for 
h{t): 



HO 



dlnV(e) 
3 V dln(l + |) 



where V(£) is given by the fitting function 
1 + 0.0158 C 2 + 0.000115 



1 + 0.926 - 0.0743 £ 3 + 0.0156 £ 4 



1/3 



(38) 



(39) 



Figure (||) shows the simulation data represented by the fit 
(solid line) ( [Hamilton et al. 1991 ) and the best fit (dashed 
line), obtained in our model, for q ~ 0.02. We note that 
the fit is better than 5% for all values of density contrast 
8 > 15. The change in the fit is very marginal if one imposes 
the boundary condition h{8) — * 1 for 8^1, instead of con- 
straining the curves to match at their peaks (for example, 




200 250 300 
Density contrast 

Figure 4. The figure shows the best fit curve for the h function 
(dashed line) to the simulation data (solid line). The simulation 
results are obtained from Hamilton et al. (1991) and the fit is 
obtained by adjusting the value of q parameter until the curves 
coincide. 



the change in the peak height is 
above condition at 8 — 10000). 



1%, if we impose the 



Figure (a) shows the plot of scaled radius y q (x) vs x, ob- 
tained by integrating equation ©, with q = 0.02. The fig- 
ure also shows an accurate fit (dashed line) to this solution 
of the form 

. . x + ax 3 + bx 5 . . 

VM = T-— ( 4 °) 



1 + cx z + bx 5 
with a — —3.6, b — 53 and c = 



-12. This fit, along with 



values for 



and 



completely specifies our model 



through equation (pl|). It can be observed that (r V i r /R„ 
is approximately 0.65. It is interesting to note that the value 
obtained for the (r„; r / Rmax) ratio is not very widely off the 
usual value of 0.5 used in the standard spherical collapse 
model, in spite of the fact that no constraint was imposed 
on this value, ab initio, in arriving at this result. Part of 
this deviation may also originate in the fit which has been 
used for h(£); Hamilton et al. (1991) noticed that objects 
virialised at Rmax/r V i r ~ 1.8, instead of 2, in their simula- 
tions. 

Finally, figure ([]) compares the non-linear density contrast 
in the modified SCM (dashed line) with that in the stan- 
dard SCM (solid line), by plotting both against the lin- 
early extrapolated density contrast, 8l- It can be seen (for 
a given system with the same z m ax and r„ir) that, at the 
epoch where the standard SCM model has a singular be- 
haviour (Sl ~ 1.686), our model has a smooth behaviour 
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1.6 
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1.2 



0.8 
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0.5 1 1.5 2 2.5 3 3.5 

x 



Figure 5. The figure shows a plot of the scaled radius of the 
shell y q as a function of scaled time x (solid line) and the fitting 
formula y q = (x + ax 3 + bx 5 )/(l + ex 3 + 6x 5 ), with a = —3.6, 
b = 53 and c = — 12 (dashed line) (See text for discussion) 




1.2 1.3 1.4 1.5 1.6 1.7 

Linear density contrast, (5 L 



Figure 6. The figure shows the non-linear density contrast in the 
SCM (solid line) and in the modified SCM (dashed line), plotted 
against the linearly extrapolated density contrast 8^ (discussion 
in text). 



with 8 w 110 (the value is not very sensitive to the exact 
value of q). This is not widely off from the value usually 
obtained from the ad hoc procedure applied in the stan- 
dard spherical collapse model. In a way, this explains the 
unreasonable effectiveness of standard SCM in the study of 
non-linear clustering. 

As mentioned earlier, deviations from spherical symmetry 
are expected to be small at early epochs and to grow as the 
system evolves. One would thus expect the two curves of fig- 
ure (^) to approach each other as 8 — > 1 (from above). Fur- 
ther, the curves should overlay in the linear regime (Sl <C 1). 
It can be seen from the figure that the 2 curves do approach 
each other as Sl reduces towards unity. However, the MSCM 
has been obtained using a Taylor expansion in (1/5); it is 
clearly not applicable for iS C 1. Further, the region 8 > 15 
has been used to fit the function h(8) to the data of Hamil- 
ton et al. (1991). Hence, one cannot compare the curves in 
the linear regime. 

Figure ([]) also shows a comparison between the standard 
SCM and the MSCM in terms of 8 values in the MSCM 
at two important epochs, indicated by vertical arrows, (i) 
When R = Rmax/Z in the SCM, i.e. the epoch at which the 
SCM virializes, <5(MSCM) ~ 83. (ii) When the SCM hits the 
singularity, (5i ~ 1.6865), <5(MSCM) ~ 110. 

We note, finally, that figure ([]), which shows the effects of 
evolution as a mapping from linear to non-linear density con- 
trasts, contains a subtle implicit assumption regarding the 
definition of the non- linear density contrast. The radius R of 
a system is not a rigorously defined quantity in the absence 
of spherical symmetry, and obviously, any argument involv- 
ing 'virialization' precludes strict spherical symmetry. It is, 
however, a conventional practice to define the 'radius', R, 
even for a virialized system without strict spherical symme- 
try. For example, this approach is used to define the density 
contrast at 'virialization' (which has the value, 5 v i r ~ 200) 
in the standard SCM. In our model we have an explicit equa- 
tion for R; once R and M are given, the non-linear density 
contrast is a well-defined quantity. 



5 RESULTS AND CONCLUSIONS 

In this paper, we have shown how the Taylor expansion of 
a term in the equation for the evolution of the density con- 
trast, 8, in inverse powers of 8, allows us to have a more 
realistic picture of spherical collapse, which is free from ar- 
bitrary abrupt "virialization" arguments. Beginning from a 
well motivated ansatz for the dependence of the "virializa- 
tion" term on the density contrast we have shown that a 
spherical collapse model will gracefully turn around and col- 
lapse to a constant radius with 8 ~ 110 at the same epoch 
when the standard model reaches a singularity. Figure (|E]) 
shows clearly that the singularity is avoided in our model due 
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to the enhancement of deviations from spherical symmetry, 
and consequent generation of strong non-radial motions. 
We derive an approximate functional form for the virializa- 
tion term starting from the physically reasonable assump- 
tion that the system reaches a constant radius. This assump- 
tion allows us to derive an asymptotic form for the virializa- 
tion term, with the residual part adequately expressed by 
keeping only the first and second order terms in a Taylor 
series in (1/S). It is shown that there exists a scaling rela- 
tion between the coefficients of the first and second order 
terms, essentially reducing the virialization term to a one 
parameter family of models. 

The form of the h function published in the literature, along 
with a tentative mapping from 8 to f , in the non-linear and 
quasi-linear regimes, allow us to further constrain our model, 
bringing it in concordance with the available simulation re- 
sults. Further, it is shown that this form for the virialization 
term is sufficient to model the turnaround behaviour of the 
spherical shell and leads to a reasonable numerical value for 
density contrast at collapse. 

There are several new avenues suggested by this work which 
we plan to pursue in the future. 

(i) The assumption hsc — > 1, -R — > r V i r is equivalent to "sta- 
ble clustering", in terms of the statistical behaviour. Since 
stable clustering has so far not been proved conclusively in 
simulations and is often questioned, it would be interesting 
to see the effect of changing this constraint to ft — » constant 
for t — > oo. 

(ii) The technique of Taylor series expansion in (1/5) seems 
to hold promise. It would be interesting to try such an at- 
tempt with the original fluid equations and (possibly) with 
more general descriptions. 

(iii) It must be stressed that we used the 8 — £ mapping 
ibly the weakest part of our analysis , conceptually 



pos.' 

nnlyl tn fW a -iralnp nf q Wp fnnlrl Viaw iigprl gnmg ViigVi 



resolution simulations to actually study the evolution of a 
realistic overdense region. We conjecture that such an anal- 
ysis will give results in conformity with those obtained here. 
We plan to investigate this — thus eliminating our reliance 
on the 8 — £ mapping — in a future work, 
(iv) Finally, the curves of figure (6) can be used to describe 
the spatial distribution of virialised haloes (see, for example, 
Mo & White 1996; Sheth 1998). It would be interesting to 
investigate how things change when the MSCM is used in 
place of the standard spherical collapse model. 



dh h 



4ttF 



9 / <91n F\ dhl 9 



where we have used the ansatz, h = h(£) (Hamilton et al. 
1991; Nityananda & Padmanabhan 1995). In the above, we 
have defined 



£(x,a) — — J dx( ! (x,a)x 



F=|| + 3(1 + C) , X=lnx 



n 



(42) 



(43) 



(44) 



where II and E are parallel and perpendicular peculiar ve- 
locity dispersions (Ruamsuwan & Fry 1992; Kanekar 1999). 
Finally, we have assumed that the 3-point correlation func- 
tion has the hierarchical form (Davis & Peebles 1977; Ru- 
amsuwan & Fry 1992) 



Cl23 = Q(Cl2^13 + £13^23 + £12^23) 

and defined 



M = / dz 



, cos 6 

£( z - x )— — 



(45) 



(46) 



In the non-linear regime, £ S> 1, the stable clusteri ng ansatz 
yields a scale-invariant power- law behaviour for £ (Davis & 



Peebles 1977), with £ tx a^-^x' 1 , if h -> 1 as £ -» 00. In 



this limit, we have 

F= (3-7)? + 3 
and 

<91n F _ 
dX 



^4: 

-7 <e 



(47) 



(48) 



Further, we can write ( Yano fc Gouda 1997 ) M = M'x£ , 
where M' is a constant. Thus, equation (fill) reduces, in the 
non-linear regime, to 



6 APPENDIX 

The zeroth and first moments of the 2 nd BBGKY equation 
(Ruamsuwan & Fry (1992), equations (22) and (23)) can be 
combined to obtain the follo wing equation for the dimen- 
sionless function h — —v/dx (Kanekar 1999) 



3/^ 



h +- 



(4 - 7) + 



3-7 
3 7 



(3 - 7)£ 



3M' 




47T 




dh\ 


h\ + 




dX 



+ 1 



• 2ft1 + O 



(49) 



where we have retained terms upto order constant in £. We 
now assume that h% and h\ are functions of £ alone, to first 
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order. This yields 



3/i£ — - h + - 
d£ 2 



with 



3M' 
4% 



+ 1 



G(0 + O(|) (50) 



G(0 



(4 - 7) + 



37 



(3-7)€ 



„ _ -dhl 



2h\C0 (51) 



Clearly, if h — > 1 as £ — > oo, we must have 



G(0 = 



3M' f- 

4tt ^ 



+ 1 



-5 + (|) .■">>> 



i.e. G(£) » -9M'£/4tt(3 - 7) for £ > 1. 
Since G(£) tends to the above asymptote at late times, the 
residual part can be expanded in a Taylor series in l/£. Re- 
taining the first two terms of the expansion in equation (|50|), 
we obtain 
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3ft£ A a + + = + ) (53) 

d£ 2 2 £ £ 

This is exactly the same as equation (|2^), with £ replacing 
5. G(£) thus plays the same role as S(8) in the stabilising 
of the system against collapse. This clearly implies that the 
velocity dispersion terms, /i| and h\, will contribute to the 
support term; we are hence justified in writing the virializa- 
tion term in the more general form 

S = a (a 2 -2Q 2 ) + f(a,x) (54) 

where f(a, x) contains contributions from effects arising 
from shell crossing, multi-streaming, etc. 
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